Kuiper Belt–like Hot and Cold Populations of Planetesimal Inclinations in the β Pictoris Belt Revealed by ALMA

Figure 1 shows the continuum emission map of the β Pictoris system obtained from the combined, naturally weighted, visibility data set. A peak signal-to-noise ratio (S/N) of 16 is reached at the location of the central star, for the rms noise level of 12 μJy/beam measured in an emission-free region of the continuum map. Emission from the belt is also clearly detected with a peak S/N of ∼14, as measured to the southwest (SW) of the star along the midplane. The spatially integrated 1.33 mm flux measured on the primary-beam-corrected map within a 16'' × 4'' box centered on the star is 20 ± 2 mJy (where the uncertainty is dominated by the 10% contribution from flux calibration; Fomalont et al. 2014).

Zoom In Zoom Out Reset image size

The stellar contribution as measured from the peak flux at its location is 190 μJy, which should strictly be considered an upper limit due to extra disk emission lying along the line of sight to the star. In order to separate the stellar flux from the disk emission, we extract the subset of the data set containing visibilities at u − v distances >300kλ. This u − v cutoff was chosen through inspection of CLEAN images, as an optimal trade-off to both ensure that all large-scale emission from the disk is filtered out, while retaining as much data as possible for maximum sensitivity. The obtained CLEAN image has an rms noise level of 18 μJy/beam for a beam size of 022 × 019, and a compact point source is detected at the stellar position. A simple 2D Gaussian fit within CASA yields a stellar flux of 80 ± 20 μJy at 1.33 mm. This detection is consistent with the expected stellar flux of ∼86 μJy obtained by extrapolating the Rayleigh–Jeans tail of the stellar blackbody-like emission from near-IR wavelengths. We therefore assume a stellar flux of 86 μJy at 1.33 mm throughout the rest of this paper.

In Figure 2, we show the profile of the continuum emission vertically integrated within ±20 au of the disk midplane, defined here as the position angle of the main disk observed in scattered light (293; Lagrange et al. 2012). Contrary to what is seen for CO emission in the same data set (green line), the millimeter-sized dust traced by the 1.33 mm continuum observations (blue line) shows no significant peak brightness asymmetry between the SW and northeast (NE) sides of the disk. This is in contrast with the SW/NE brightness asymmetry, between 30 and 80 au on either side the star, reported at 885 μm by Dent et al. (2014). The SW/NE ratio measured within the same radial distances as Dent et al. (2014), vertically integrated between ±20 au, is 1.03 ± 0.04 in the new 1.33 mm data set, consistent with no asymmetry. The same SW/NE ratio is 1.05 ± 0.03 in the 885 μm archival data set, indicating again a negligible (5 ± 3)% asymmetry in this radial region. As this measurement can be affected by imaging in the presence of incomplete u − v sampling and subsequent deconvolution via the CLEAN algorithm, we postpone further discussion of this issue in the context of visibility models to Section 4.6.

Zoom In Zoom Out Reset image size

Figure 3 (top) shows vertical profiles of emission perpendicular to the disk midplane as a function of distance to the star, which were obtained by integrating emission every 15 au along the radial (midplane) direction. As shown by comparing the profiles (blue shaded region) to the resolution element of our observation (black dotted line), the dust disk is clearly resolved in the vertical direction. Following the procedure described in detail in Section 3.2 and Appendix B of Matrà et al. (2017), we fit Gaussians (purple dashed lines) to the vertical profiles at each midplane location, and derive vertical best-fit Gaussian centroids (forming the disk spine, Figure 3 middle panel) and widths (measured as FWHM, Figure 3 bottom panel). We report two main findings as follows.

• 1.  
  In contrast with the CO J = 2–1 line emission, we find that dust emission is not well represented by a single vertical Gaussian, being more centrally peaked and having broader wings extending out to ±40 au from the sky-projected midplane. This is evident in Figure 4, which shows the residual vertical profile obtained by integrating emission across the entire midplane after vertically centering it so that the centroid of the emission aligns with the midplane. The best-fit Gaussian profile (red, with a standard deviation of 8.2 ± 0.2 au) leaves significant residuals (Figure 4, bottom), whereas a double-Gaussian profile (green, with standard deviations of 5.1 ± 0.3 and 15.7 ± 1.4 au) considerably improves the fit, leaving no significant residual emission. We find that relative flux calibration systematics of ±10% between the data sets with compact and extended baselines (tracing vertically broader and narrower emission, respectively) cannot account for this discrepancy, which must therefore be physical in origin. This shows that the double-Gaussian or Lorentzian profiles needed to fit scattered light observations (Golimowski et al. 2006; Lagrange et al. 2012; Apai et al. 2015) reflect the vertical structure of the millimeter emission that traces the dust-producing planetesimals.
• 2.  
  The belt spine (Figure 3, middle) presents little vertical displacement along the midplane and is largely consistent with the PA of the main disk observed in scattered light. This is significantly different from the CO J = 2–1 emission, whose spine presents an extra tilt angle dPA of ∼4° similar to that of the warp interior to ∼80 au seen in scattered light (Apai et al. 2015; Matrà et al. 2017). Therefore, we conclude that an inner disk tilt as large as observed in scattered light and CO observations is not detected in the structure of millimeter grains. Assuming the mm grains are tracing the planetesimal distribution, this shows that planetesimals are not perturbed into a warp—or at least not into one as pronounced as that observed in the distribution of the smallest grains observed in scattered light. At the same time, the tilt present in both the scattered light and the CO distribution, with the CO clump clearly offset by 4–6 au above the main disk midplane, suggests a connection between the two.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We model the β Pic belt using mm dust density distributions described in the following subsections, and the star as a point source with a flux density of 86 μJy located in the geometrical center of the circular belt. For each model realization, we solve the radiative transfer at 1.33 mm using RADMC-3D,⁸ assuming the dust radial temperature profile expected from a blackbody around a star of 8.7 L_⊙ (Kennedy & Wyatt 2014). Since we are interested in the structure rather than the overall amount of material in the belt, we fix the dust opacity and fit for the belt mass as a normalization factor, noting that this does not affect our conclusions on the belt structure. In addition to the parameters that govern the belt dust density distribution and its overall mass, we fit for belt geometry parameters (inclination I and position angle PA).

Having created the model belt image, we produce visibility data separately for each of the ALMA arrays (12 m and ACA), for each configuration of the 12 m array and for each of the two mosaic pointings used for the compact configuration observations. We do so by first shifting the model belt image by R.A. and decl. offsets that we leave as free parameters in the fit, and multiplying it by the primary-beam response of the relevant array, at its appropriate sky location for each observation. We then Fourier transform the image to produce a grid of model complex visibilities, which we evaluate at the u − v locations sampled by the observation in question as described, for example, in Marino et al. (2016) and Tazzari et al. (2018). Finally, we add a point source representing the star to the model visibilities, located in the geometric center of the belt, and with its flux appropriately attenuated, given its location with respect to the primary-beam center.

The model visibilities are then fitted to the observed visibilities through a Monte Carlo Markov Chain (MCMC) method implemented through the emcee package (Foreman-Mackey et al. 2013). Eight nuisance parameters for the spatial shifts (in R.A. and decl., so two for each data set), two parameters to determine the viewing geometry (I and PA), one normalization parameter setting the mass of the belt for a fixed opacity, and n parameters governing the dust distribution make up the 11 + n-dimensional parameter space that we explore in our fit. Through the emcee package, we sample the 11 + n-dimensional posterior probability distribution of the explored parameters using the affine-invariant MCMC ensemble sampler of Goodman & Weare (2010). The posterior probability distribution is computed assuming uniform priors on all of the probed parameters, and a likelihood function proportional to ${e}^{-{\chi }^{2}/2}$, where χ² is the usual definition of the chi-square function. The MCMC was run for each tested model with a number of walkers equal to 10× the number of free parameters and more than 2500 steps, ensuring each of the chains had reached convergence.

The uncertainty σ on each visibility data point is defined using the visibility weights w delivered by the ALMA observatory, where ${\sigma }_{{u}_{i}-{v}_{i}}=1/\sqrt{{w}_{{u}_{i}-{v}_{i}}}$. While assuming the delivered weights and hence the uncertainties are correct for different u − v points relative to each other, we left a constant rescaling factor, equal for all u − v points within a given data set, as a free parameter in our model fit of Section 4.2. This is justified by previous ALMA modeling studies also needing to rescale the weights to appropriately represent the true scatter in the visibilities (e.g., Marino et al. 2018). We find this rescaling factor to be different but very well constrained to 0.4508 ± 0.0006, 0.927 ± 0.014, 0.690 ± 0.003, and 0.698 ± 0.003 for the 12 m extended, ACA, and for the two pointings of the 12 m compact configuration, respectively. We therefore fix the rescaling factors to these values before producing the image in Figure 1 and do the same for all subsequent model fitting runs.

We begin by assuming the simplest dust density distribution model with the minimum number of free parameters—namely a radially and vertically Gaussian, axisymmetric belt with radially constant aspect ratio. The distribution of the dust density ρ in the belt reads

$$\begin{eqnarray}&&\rho ={\rho }_{0}\ {e}^{-\tfrac{{(r-{r}_{{\rm{c}}})}^{2}}{2{\sigma }^{2}}}\displaystyle \frac{{e}^{-\tfrac{{z}^{2}}{2{({hr})}^{2}}}}{\sqrt{2\pi }{hr}},\end{eqnarray} \tag{ 1 }$$

where r and z are cylindrical coordinates, r_c and σ are the center and standard deviation of the radially Gaussian belt, h is the aspect ratio (here independent of radius), and ρ₀ is a normalization factor proportional to the total dust mass in the belt, which we fit for.

The best-fit parameters of this model, defined as the ${50}_{-34}^{+34}$th percentiles of their marginalized posterior probability distributions, are listed in Table 1 (Model 1). Within the framework of this model, and in qualitative agreement with Figure 1, we find the disk to be near edge-on, with an inclination from face-on of ${86.6}_{-0.3}^{+0.4}$ degrees, and with its Gaussian radial distribution being centered at ∼105 au with a width of ∼92 au.

Table 1.  Best-fit Model Parameters

Parameter	Unit	Model 1	Model 2	Model 3	Model 4	Model 5
rc	au	${104.9}_{-1.1}^{+1.1}$	${107.5}_{-1.0}^{+1.0}$	${106.3}_{-1.3}^{+1.4}$	${107.8}_{-1.0}^{+1.0}$	${105.2}_{-1.4}^{+1.4}$
σ	au	${38.9}_{-1.1}^{+1.3}$	${36.7}_{-1.3}^{+1.0}$	${36.6}_{-1.2}^{+1.1}$	${36.4}_{-1.3}^{+1.1}$	${36.3}_{-1.2}^{+1.1}$
h		${0.070}_{-0.004}^{+0.004}$	${0.014}_{-0.006}^{+0.006}$	${}^{* }{0.020}_{-0.008}^{+0.010}$	${0.054}_{-0.004}^{+0.004}$	${}^{* }{0.072}_{-0.012}^{+0.014}$
I	◦	${86.6}_{-0.3}^{+0.4}$	${88.6}_{-0.3}^{+0.2}$	${88.8}_{-0.3}^{+0.4}$	${89.3}_{-0.5}^{+0.4}$	${89.5}_{-0.4}^{+0.3}$
PA	◦	${30.0}_{-0.1}^{+0.1}$	${29.7}_{-0.1}^{+0.1}$	${29.7}_{-0.1}^{+0.1}$	${29.8}_{-0.1}^{+0.1}$	${29.7}_{-0.1}^{+0.1}$
h2		⋯	${0.110}_{-0.006}^{+0.008}$	${}^{* }{0.14}_{-0.02}^{+0.02}$	⋯	⋯
a		⋯	${0.20}_{-0.04}^{+0.05}$	${0.21}_{-0.04}^{+0.05}$	⋯	⋯
β		⋯	⋯	${0.7}_{-0.2}^{+0.2}$	⋯	${0.4}_{-0.2}^{+0.2}$
p		⋯	⋯	${0.7}_{-0.2}^{+0.2}$	${0.9}_{-0.1}^{+0.1}$	${0.85}_{-0.07}^{+0.07}$
Free Parameters		14	16	17	15	16
${\chi }^{2}$		1456364.9	1456297.6	1456297.4	1456312.2	1456297.3
Δχ2		⋯	−67.3a	−0.2b	+14.6b	−0.3b
ΔAIC		⋯	−63.3a	+1.8b	+12.6b	−0.3b
ΔBIC		⋯	−38.9a	+14.0b	+0.33b	−0.3b

Notes. Model 1: radially Gaussian ring with peak r_c and standard deviation σ, with a single-component, Gaussian vertical profile with a radially constant aspect ratio h. Model 2: same as model 1, with a second, broader Gaussian component to the vertical profile, with aspect ratio h₂. The ratio of the peak value of the second component to that of the first component is a. Model 3: same as model 2, but with a radially varying scale height for both Gaussian components, following H ∝ r^β, or equivalently h ∝ r^β−1. Model 4: same as model 1, but with a vertical profile described by a non-Gaussian functional form (where $\rho \propto \exp \left[-{\left(\tfrac{| z| }{\sqrt{2}{hr}}\right)}^{p}\right]$, see Equation (3)), with a radially constant aspect ratio h. Model 5: same as model 4, but with a radially varying aspect ratio h ∝ r^β−1. For Model 3 and 5, * symbols indicate that aspect ratios are quoted at a radius of 50 au.

^aDifference with respect to Model 1. ^bDifference with respect to Model 2.

Download table as:  ASCIITypeset image

Figure 5, top left, shows the residuals obtained after subtracting the best-fit model visibilities from the data and imaging the residual visibility data set. The bottom left panel of Figure 5 shows the residual radial profile constructed as in Figure 2; the lack of significant residuals and SW/NE asymmetries reaffirms our conclusion that no clear departure from axisymmetry is present, within our measurement uncertainties, in the ALMA 1.33 mm data set.

Zoom In Zoom Out Reset image size

On the other hand, the right panel shows the vertical profile obtained after radially integrating the belt's residual emission as done in Figure 4. This displays significant residuals very similar to those obtained by simply subtracting a single Gaussian from the vertical flux distribution in the image plane. This confirms that the vertical distribution of mm dust in the planetesimal belt around β Pictoris cannot be reproduced by a single Gaussian.

Motivated by the appearance of the vertical emission distribution in the image plane, we construct a belt model with a vertical distribution characterized by a narrower and broader Gaussian. The Gaussians are both centered in the belt's midplane (z = 0) and are characterized, respectively, by aspect ratios h and h₂ leading to standard deviations (scale heights) hr and h₂r at any given radius r. They comprise a fraction a and 1 − a of the total belt mass. The functional form of the dust density distribution thus reads

$$\begin{eqnarray}&&\rho ={\rho }_{0}\ {e}^{-\tfrac{{(r-{r}_{{\rm{c}}})}^{2}}{2{\sigma }^{2}}}\left[\displaystyle \frac{a}{\sqrt{2\pi }{hr}}{e}^{-\tfrac{{z}^{2}}{2{({hr})}^{2}}}+\displaystyle \frac{1-a}{\sqrt{2\pi }{h}_{2}r}{e}^{-\tfrac{{z}^{2}}{2{({h}_{2}r)}^{2}}}\right].\end{eqnarray} \tag{ 2 }$$

This model is significantly better at describing the vertical distribution of the belt, as seen when comparing its residuals (Figure 6) to the single vertical Gaussian model (Figure 5). Comparing best-fit models (Model 2 and 1, respectively, in Table 1) through an Akaike or Bayesian Information Criterion (AIC, or BIC), which penalize models with a larger number of free parameters, clearly favors the model with the more complex vertical structure. Formally, the large BIC differential (or analogously the ratio of the model likelihoods) indicates the model with a double-Gaussian vertical structure is >10⁸ times more likely to be better at describing the data than the model with a single Gaussian, providing strong evidence in support of our conclusion.

Zoom In Zoom Out Reset image size

The narrow Gaussian component has a vertical FWHM of ∼2 au at 150 au, which is similar to the size of our resolution element and only resolved through the forward modeling presented here. Note also that the inclination of the double-Gaussian model is significantly higher (closer to edge-on) than found for the single-Gaussian model. This is likely attributable to the single-Gaussian model being driven to a less edge-on configuration in an attempt to capture the broad sky-projected vertical distribution of the disk emission. The position angle east of north of the mm dust disk (${29.7}_{-0.1}^{+0.1}$) is marginally consistent within the uncertainties with the outer disk observed in scattered light (${29.3}_{-0.1}^{+0.2}$; Lagrange et al. 2012).

The appearance of the residuals in Figure 6, after subtracting a model where the aspect ratio is the same at all belt radii, indicates that we have already achieved a good fit to the data. Nonetheless, we test whether we can detect a radial dependence of the vertical aspect ratio through a model where the scale height is a function of radius, and this dependence is parameterized through a flaring parameter β. The dust density distribution reads the same as in Equation (2), but with h and h₂ being a function of radius following $h(r)={h}_{50\mathrm{au}}{\left(\tfrac{r}{50\mathrm{au}}\right)}^{\beta -1}$ and similarly ${h}_{2}(r)={h}_{\mathrm{2,50}\mathrm{au}}{\left(\tfrac{r}{50\mathrm{au}}\right)}^{\beta -1}$. Then, β = 1 indicates a constant aspect ratio h = H(r)/r as considered in the previous models, β > 1 a flared disk, and β = 0 a disk with a radially constant scale height.

The marginalized posterior probability distribution of β results in $\beta ={0.7}_{-0.2}^{+0.2}$, indicating that models with a scale height increasing with radius (β > 0) are significantly preferred, with a marginal improvement achieved for models with an aspect ratio decreasing with radius (β < 1). However, we find no evidence that the extra free parameter β leads to an overall improvement in the fit to the data, which is confirmed by visual inspection of the residuals. We therefore conclude that our data supports models where the scale height increases with radius, with a flaring parameter β consistent with one.

Our results do not rule out that a different parametric prescription for the vertical density distribution could describe the data as well as the double-Gaussian parameterization adopted here. For example, we here attempt to use a model where the drop-off of the vertical density distribution differs from a Gaussian (as implemented by Ahmic et al. 2009). The functional form reads

$$\begin{eqnarray}&&\rho ={\rho }_{0}\ {e}^{-\tfrac{{(r-{r}_{{\rm{c}}})}^{2}}{2{\sigma }^{2}}}{{Ce}}^{-{\left(\tfrac{| z| }{\sqrt{2}{hr}}\right)}^{p}},\end{eqnarray} \tag{ 3 }$$

with $h(r)={h}_{50\mathrm{au}}{\left(\tfrac{r}{50\mathrm{au}}\right)}^{\beta -1}$ (as in Section 4.4), and $C\,={\left[{\int }_{-\infty }^{+\infty }{e}^{-{\left(\tfrac{| z| }{\sqrt{2}{hr}}\right)}^{p}}\right]}^{-1}$ ensuring normalization of the vertical distribution. Distributions with p < 2 allow the vertical distribution to drop faster than a Gaussian the further we move from the midplane, but to have broader wings; Ahmic et al. (2009) found this scenario to be preferred in scattered light observations of the β Pic belt, with p values around or below 1.

We tried two model runs with the above functional form and p left as a free parameter. In the first (Model 4), we fix the aspect ratio (β = 1) as done for the Model 1 and 2 runs. In the second run (Model 5), we also leave the flaring parameter β free. Our fitting results for both model runs, shown in Table 1, indicate that p ∼ 0.9, in broad agreement with the scattered light findings, and consistent with the broad-winged distribution derived from the image analysis.

We also find that best-fit Model 4 and 5 (with p = 0.9) are significantly better than a single Gaussian (Model 1). Compared with the double Gaussian, constant aspect ratio model (Model 2, see Table 1, bottom), we find that Model 4 (with β fixed to 1) is slightly worse, whereas Model 5 (with β ∼ 0.4) is just as good at reproducing the data, with extremely similar residuals. In summary, there is no evidence in the mm data that a model with a single population of mm grains that has a vertical drop-off rate that is shallower than Gaussian (Equation (3) with p ∼ 0.9, Models 4 and 5) is preferred to a double-Gaussian model (Model 2). In fact, the convergence of Models 4 and 5 to a broad-winged distribution confirms our main finding that there exists a population of high inclination (dynamically hot) particles in the β Pic disk. For later interpretation, we prefer to use the double Gaussian simply to facilitate direct comparison with the Kuiper Belt (see Section 5.2.3).

Independent of the models explored previously, the vertically integrated radial profile of the residuals of the 1.33 mm data shows no significant SW/NE asymmetry (bottom left panel of Figures 5 and 6, where the shaded region is the ±1σ confidence interval). Averaging emission between 30 and 80 au on each side of the star, as in Section 3, we find a negligible residual flux difference between the NE and the SW of 82 ± 85 μJy. Varying the ranges over which we radially average does not affect this result. We conclude that we do not detect a significant SW/NE asymmetry in the new 1.33 mm data set.

Given the discrepancy with the 15% asymmetry reported in the 885 μm data set (Dent et al. 2014), we independently re-reduced the data set as described in Section 2.2 of Matrà et al. (2017). Just like the 1.33 mm, compact configuration data set, this consisted of a mosaic pointed ±5'' from the stellar location and along the belt's PA. We imaged the continuum using the tclean task (CASA v5.1.0) in multi-frequency synthesis, mosaic mode, using natural weights and multiscale deconvolution to best recover faint, extended emission. We obtained an image with synthesized beam size of 070 × 055 at a PA of 644, and measure an rms noise level of 70 μJy beam⁻¹. The spatially integrated line flux measured on the primary-beam-corrected map as in the 1.33 mm data, is 66 ± 6 mJy.

As done for the 1.33 mm data, we extract the subset of the data set containing visibilities at u − v distances >300kλ to separate the stellar emission from the disk contribution. We find that the star is not significantly detected. This allows us to set a 3σ upper limit of <786 μJy on the stellar flux at 885 μm, which is consistent with the expected 196 μJy obtained through Rayleigh–Jeans extrapolation from near-IR wavelengths. We therefore include the star with a fixed flux of 196 μJy in the modeling that follows.

We simultaneously fit the visibilities obtained for both pointings as described in Section 4.1, using the single and double vertical Gaussian models described in Sections 4.2 and 4.3. For each of the two models, we fix the model parameters to the best-fit values found for the higher resolution and sensitivity, 1.33 mm data set (Table 1, Model 1 and 2). We leave as free parameters the scale factor ρ₀ governing the total flux of the belt, which is different at the shorter wavelength, as well as an R.A. and decl. offset and a weight-rescaling factor for each of the two pointings that should be different from the 1.33 mm data sets.

Figures 7 and 8 show residual maps for the single and double vertical Gaussian models, respectively. We find the residuals, particularly the radial and vertical profiles, to be remarkably similar to those obtained for the 1.33 mm observations. In fact, despite the lower resolution and sensitivity, the 885 μm data set is also better fitted by a model with a double over a single-Gaussian vertical distribution, and also shows no significant SW–NE asymmetry. Averaging emission between 30 and 80 au on each side of the star, as in Dent et al. (2014), we find a negligible residual flux difference between the NE and the SW sides of 60 ± 130 μJy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An interesting hypothesis is that the two Gaussian populations have different spectral indices, which could be caused by a difference in the size distribution of grains between the two populations. While this can in principle be tested with resolved observations at different wavelengths, such as presented here, the coarse resolution of the 885 μm data (which, at 100 au from the star, is of order the scale height of the broad Gaussian population) and the relatively small difference in wavelength between the two data sets do not allow us to test this hypothesis. Future ALMA observations at a S/N and resolution comparable to the 1.33 mm data set presented here, and at a wavelength as far as possible from 1.33 mm, are needed to explore this idea.

The observed inclination distributions of planetesimals in the Kuiper Belt was first investigated by Brown (2001, B01 hereafter). The authors used a Gaussian multiplied by a sine function to model the observed inclinations of Kuiper Belt Objects, and showed that this functional form naturally arises from N-body simulations from gravitational perturbations, or stirring, within a planetesimal disk with initially zero inclinations. The function reads $f(i){di}\propto \sin i\ {e}^{-\tfrac{{i}^{2}}{2{\sigma }_{i}^{2}}}{di}$ and for small angles corresponds to a Rayleigh distribution, $f(i)\propto \tfrac{i}{{\sigma }_{i}^{2}}{e}^{-\tfrac{{i}^{2}}{2{\sigma }_{i}^{2}}}$. This result is consistent with other dynamical studies also finding that gravitational perturbations of planetesimals in a thin disk lead to Rayleigh distributions of inclinations and eccentricities (e.g., Ida & Makino 1992; Lissauer & Stewart 1993). For direct comparison with the Kuiper Belt, we here follow the notation of B01 but note that dynamical studies typically express the inclination distribution as a function of the mean of the squares of the inclinations $\langle {i}^{2}\rangle$, with $\langle {i}^{2}\rangle =2{\sigma }_{i}^{2}$ in the expressions quoted previously. Figure 9 (left) shows an example inclination distribution (following B01) for a population of particles with a mean inclination of ∼79 (orange points) and shows how this is well matched to a Rayleigh distribution (blue line).

Zoom In Zoom Out Reset image size

Under the assumption of circular orbits, this inclination distribution of planetesimal orbits f(i) in a belt at a given radius can be linked to the distribution of latitudes L(θ) (about a planetary system's ecliptic plane) at which the planetesimals are found (B01), through

$$\begin{eqnarray}&&L(\theta )={\int }_{\theta }^{\pi /2}f(i)\displaystyle \frac{\cos \theta }{\sqrt{{\sin }^{2}i-{\sin }^{2}\theta }}{di}.\end{eqnarray} \tag{ 4 }$$

The integral can be simplified through use of the small angle approximation (sinx ∼ x for both i and θ), given the small inclination angles considered here. For example, Figure 9 (center and right) shows the latitudinal and height distribution derived from an inclination distribution (left) with a mean $\bar{i}$ of 79 and an rms of 89, corresponding to the hot population of the β Pic belt. At the height corresponding to 3σ of the distribution (16.5 au, implying a latitude of 189), the small angle approximation is accurate to within 2%. Taking the small angle approximation, the integral in Equation (4) becomes

$$\begin{eqnarray}&&L(\theta )\sim \displaystyle \frac{\cos \theta }{{\sigma }_{i}^{2}}{\int }_{\theta }^{\pi /2}\displaystyle \frac{{{ie}}^{-\tfrac{{i}^{2}}{2{\sigma }_{i}^{2}}}}{\sqrt{{i}^{2}-{\theta }^{2}}}{di},\end{eqnarray} \tag{ 5 }$$

and can be solved analytically leading to the expression

$$\begin{eqnarray}&&L(\theta )=\sqrt{\displaystyle \frac{\pi }{2{\sigma }_{i}^{2}}}\cos \theta \ \mathrm{erf}\left(\sqrt{\displaystyle \frac{\tfrac{{\pi }^{2}}{4}-{\theta }^{2}}{2{\sigma }_{i}^{2}}}\right){e}^{-\tfrac{{\theta }^{2}}{2{\sigma }_{i}^{2}}}.\end{eqnarray} \tag{ 6 }$$

Given the small inclinations i expected to inevitably lead to small latitudes θ, we can again use the small angle approximation to obtain $\cos \theta \sim 1$ and $\mathrm{erf}\left(\sqrt{\tfrac{\tfrac{{\pi }^{2}}{4}-{\theta }^{2}}{2{\sigma }_{i}^{2}}}\right)\sim 1$ (the latter easily shown given, e.g., θ ≪ 5σ_i and small θ). Then, we come to the conclusion that the latitudinal distribution of particles with a Rayleigh distribution of inclinations (Figure 9, left) is simply a Gaussian with standard deviation equal to the σ_i parameter of the inclination distribution (i.e., $L(\theta )\propto {e}^{-\tfrac{{\theta }^{2}}{2{\sigma }_{i}^{2}}}$). We verify this in Figure 9 (center) by comparing an evaluation of the numerical integral in Equation (4) with the Gaussian derived analytically previously.

Finally, for small angles the distribution of heights (rather than latitudes) of the planetesimals will also be very well approximated by the Gaussian

$$\begin{eqnarray}&&\rho (z)\propto {e}^{-\tfrac{{z}^{2}}{2{\sigma }_{z}^{2}}},\end{eqnarray} \tag{ 7 }$$

where ${\sigma }_{z}=H={hr}\sim {\sigma }_{i}r$. Thus, using a Gaussian vertical distribution in our model fits of Section 4 is justified by the expectation of a Rayleigh distribution of inclinations for a belt whose dynamics is dominated by gravitational stirring from bodies within the belt alone (e.g., Lissauer & Stewart 1993; Stewart & Ida 2000). The aspect ratio of an observed vertical Gaussian distribution can then be used to derive the mean $\bar{i}$ and root mean square ${\mathrm{rms}}_{i}=\sqrt{\langle {i}^{2}\rangle }$ of the inclination distribution of belt particles, following

$$\begin{eqnarray}&&\bar{i}=\sqrt{\displaystyle \frac{\pi }{2}}h\ \mathrm{and}\ {\mathrm{rms}}_{i}=\sqrt{2}h.\end{eqnarray} \tag{ 8 }$$

In Sections 3 and 4 we showed how the vertical distribution of mm dust in the β Pictoris belt is best described by a sum of two Gaussian distributions. The derivation provided then indicates that the corresponding distribution of particle inclinations is bimodal, corresponding to the sum of two Rayleigh distributions (Figure 10). From now on, we will refer to these as the hot and cold populations of particles, with high and low inclinations, respectively. The best-fit rms inclination derived from Model 2 is ${0.156}_{-0.009}^{+0.011}$ rad (${8.9}_{-0.5}^{+0.7}$ degrees) for the hot population and ${0.02}_{-0.01}^{+0.01}$ rad (${1.1}_{-0.5}^{+0.5}$ degrees) for the cold population.

Zoom In Zoom Out Reset image size

5.2.1. External Stirring: β Pic b

The only known body that could contribute to structure in the belt is the directly imaged, super-Jupiter mass planet β Pic b. The most recent orbital determination for the planet indicates a best-fit semimajor axis of 8.8 au, an inclination to the line of sight I_p of 893 and a position angle Ω_p of 321 (Lagrange et al. 2019). When comparing to the belt as resolved by ALMA, we confirm that the orbital plane of the planet and the belt are misaligned with respect to one another. This is thanks to our high-resolution measurement of the belt's position angle Ω_b and inclination to the line of sight I_b from thermal emission, which for the first time are free of assumptions regarding the unknown phase function of the grains, as had affected previous determinations from scattered light observations.

The misalignment of the orbital plane of the planet and that of the disk can be characterized by the angle dI, the inclination of the planet orbit with respect to the disk midplane. A simple transformation between reference frames shows that

$$\begin{eqnarray}&&\cos ({dI})=\cos ({I}_{{\rm{b}}})\cos ({I}_{{\rm{p}}})+\sin ({I}_{{\rm{b}}})\sin ({I}_{{\rm{p}}})\cos ({{\rm{\Omega }}}_{{\rm{p}}}-{{\rm{\Omega }}}_{{\rm{b}}}).\end{eqnarray} \tag{ 9 }$$

We find that the belt's inclination I_b to the line of sight (see Table 1) for all the best-fit models to be consistent with that of β Pic b's orbit (I_p ∼ 893, but see Figure 5 in Lagrange et al. 2019, for an estimate of the uncertainty). If confirmed by future observations, this would indicate a planet-belt orbital plane misalignment dI ∼ 24, approximately equal to the difference in position angle between the belt and the planet's orbital plane.

Assuming this misalignment has been retained for the ∼23 Myr age of the system, the misaligned orbital plane of β Pic b should significantly affect the belt's vertical structure through secular perturbations affecting the belt's planetesimals. In particular, β Pic b should impose a forced inclination i_p = dI ∼ 24 on the belt's particles and cause their inclinations to oscillate between 0 and 2i_p (e.g., Dawson et al. 2011; Nesvold & Kuchner 2015). Particles further out in the belt are affected on longer timescales; this has long been considered the origin of the belt's warp observed in the inner regions out to ∼100 au from scattered light observations (Mouillet et al. 1997; Augereau et al. 2001), which trace the smallest grains in the collisional cascade.

In this scenario, we would expect (1) the warp to be also present for mm grains and (2) the vertical distribution of dust density, and hence mm emission, to be enhanced at 0 and 2i_p where particles spend most of their time during their oscillations. An example of the expected disk morphology imparted by β Pic b is shown in Figure 1 of Nesvold & Kuchner (2015). This shows clear vertical flux enhancements displaced ∼ ± 9 au from the midplane at 100 au from the star, which are clearly ruled out by the disk spine and vertical profiles derived from the new ALMA data (Figure 3, central and top panel). Indeed, the mm continuum emission is centrally peaked and displaced by no more than 5 au vertically (at the 3σ level) throughout the disk midplane. This would limit any sky-projected warp, if present, to less than 28 (when measured at 100 au), which roughly corresponds to the forced inclination imposed by the planet.

In summary, we find that the disk spine and centrally peaked vertical emission profile are inconsistent with models of the secular perturbations imposed by β Pic b on the belt. Taken alone, this could indicate that (1) the planet is less inclined compared to the belt than currently believed; (2) another massive planet is present in the outer regions of the belt, dominating its dynamics; or (3) β Pic b has only recently been put on a misaligned orbit and has not had time to interact with the outer belt. Such an alternative picture would have to be reconciled with the scattered light and CO inner disk tilt, and with the hot and cold inclination populations reported here.

5.2.2. Gravitational Stirring within the Belt and Constraints on the Size of the Belt's Largest Bodies

The strongest constraint on the double-Gaussian vertical distribution of the belt comes from its radial outer edge, which is least affected by line-of-sight integration of emission originating at different orbital radii. Even if the belt is misaligned with β Pic b's orbit, this outer region has not had enough time to be affected by secular perturbations from the planet within the age of the system. Therefore, any structure in this region, and particularly the hot population, is either inherited from the gas-rich protoplanetary phase of evolution (e.g., Walmswell et al. 2013) or must be the result of dynamical excitation from a large body/bodies other than β Pic b.

We here consider gravitational stirring by large bodies within the belt itself. Encounters with such bodies will act as a source of dynamical heating, setting the velocity dispersion of particles, and producing a single, Rayleigh distribution of inclinations and eccentricities (see Kokubo & Ida 2012, and references therein). Given that the observed inclination distribution of the β Pic belt clearly deviates from a single Rayleigh distribution, we can rule out gravitational stirring from within the belt alone as the source of the observed vertical structure. Nonetheless, here we consider gravitational stirring as a potential source of either the hot or cold populations, while keeping in mind that other mechanisms are required to produce the other population.

Gravitational stirring leads to a relative velocity of particles of the form ${v}_{\mathrm{rel}}=\sqrt{1.25\langle {e}^{2}\rangle +\langle {i}^{2}\rangle }{v}_{\mathrm{Kep}}$, where v_Kep is the local Keplerian velocity (e.g., Lissauer & Stewart 1993). Since $\sqrt{\langle {e}^{2}\rangle }=2\sqrt{\langle {i}^{2}\rangle }$ (Ida & Makino 1992), the rms of the inclination distribution is a direct measure of the relative velocity of particles, following

$$\begin{eqnarray}&&{v}_{\mathrm{rel}}=\sqrt{6\langle {i}^{2}\rangle }{v}_{\mathrm{Kep}}=103.9{M}_{\star }^{0.5}{{hr}}^{-0.5},\end{eqnarray} \tag{ 10 }$$

where h is unitless, r is in astronomical unit, M_⋆ is in Solar masses, and the resulting v_rel is in km s⁻¹. Then, for an assumed stellar mass of 1.75 M_⊙ (Crifo et al. 1997), the cold and hot populations in the β Pic belt (Model 2) indicate approximate relative velocities of 0.27 and 2.1 km s⁻¹ at 50 au, dropping to 0.16 and 1.23 km s⁻¹ at 150 au.

Viscous stirring from a large body (or bodies) will act to increase the velocity dispersion of planetesimals, and consequently their inclinations, over time. Following Ida & Makino (1993), assuming the system has been evolving from a perfectly cold disk (e, i ∼ 0 at t = 0 Myr) for the age of the system (t ∼ 23 Myr), we can then relate the measured inclinations to the product of the mass and the surface density of the largest bodies.

We assume we are in the dispersion-dominated regime and that stirring is dominated by large body/bodies rather than the small planetesimals themselves. Combining Equations (4.2) and (4.9) of Ida & Makino (1993) and solving for the inclination as a function of time, we obtain

$$\begin{eqnarray}&&\sqrt{\langle {i}^{2}\rangle }={\left(2{C}_{{\rm{I}}}\displaystyle \frac{M}{{M}_{\star }^{2}}{r}^{2}{\rm{\Omega }}{\rm{\Sigma }}t\right)}^{\tfrac{1}{4}},\end{eqnarray} \tag{ 11 }$$

where the constant C_I ∼ 2 (Ida & Makino 1993), Ω is the Keplerian frequency in rad s⁻¹, M_⋆ is the stellar mass in kg, t is time in s, and Σ is the surface mass density of the large stirring bodies (in kg m⁻²), each of mass M (in kg). The measured rms inclination of a planetesimal population then allows us to set a constraint on the mass times the surface density of the largest bodies of

$$\begin{eqnarray}&&M{\rm{\Sigma }}=4.43\times {10}^{3}{M}_{\star }^{\tfrac{3}{2}}{r}^{-\tfrac{1}{2}}\displaystyle \frac{{\sqrt{\langle {i}^{2}\rangle }}^{4}}{t}\ \ {M}_{\oplus }^{2}\ {\mathrm{au}}^{-2},\end{eqnarray} \tag{ 12 }$$

(where M_⋆ is now in M_⊙, r is in astronomical unit, and t is in Myr). This leads to a value of 5.9 × 10⁻⁶ ${M}_{\oplus }^{2}$ au⁻² for the cold population, and 0.022 ${M}_{\oplus }^{2}$ au⁻² for the hot population, at 150 au in the β Pic belt. This corresponds to the purple lines in the Σ versus M space shown in Figure 11.

Zoom In Zoom Out Reset image size

Producing the observed excitation however requires the bodies causing the stirring to not collide (and hence be destroyed) within the age of the system. The rate of mutual collisions for a population of bodies of equal mass M can be estimated as

$$\begin{eqnarray}&&{R}_{\mathrm{col}}=\sigma {{nv}}_{\mathrm{rel}}\left(1+\displaystyle \frac{{v}_{\mathrm{esc}}^{2}}{{v}_{\mathrm{rel}}^{2}}\right),\end{eqnarray} \tag{ 13 }$$

where σ = 4πR² is the effective cross section of the collision (with R = (0.24M/ρ)^1/3 being the radius of the large bodies, and ρ their material density), n is the number density of such bodies, and ${v}_{\mathrm{esc}}=\sqrt{{{GM}}_{\star }/R}$ is their escape velocity. As in Equation (8), we assume ${v}_{\mathrm{rel}}=\sqrt{6\langle {i}_{p}^{2}\rangle }$, where $\langle {i}_{p}^{2}\rangle$ is the rms inclination of the large bodies doing the stirring, which is likely primordial and here assumed to be 0.05. Setting the collision timescale (${t}_{{\rm{c}}}={R}_{\mathrm{col}}^{-1}$) to be greater than the age of the system, we can obtain an upper limit to the surface density of the large bodies given their mass, shown as the orange downward triangles in Figure 11. In log–log space, the slope of this surface density upper limit transitions from positive at the low mass end where v_esc ≪ v_rel, to negative at the high mass end where v_esc ≫ v_rel.

The intersection between the orange and purple lines indicates the minimum mass necessary for large bodies to produce the observed level of stirring but at the same time avoid potentially destructive collisions in 23 Myr. It can be shown analytically by combining Equations (9) and (11) at the high mass end (where v_esc ≫ v_rel) that requiring the collision timescale to be equal to the viscous stirring timescale leads to the planetesimals' relative velocities to be not exactly equal to but of order the escape velocity of the largest bodies (blue dashed–dotted line). This means that this minimum mass can be approximately estimated by simply setting the planetesimals' relative velocities to be equal to the stirring body's escape velocity.

A further requirement for the surface density of large bodies in the belt is that at least one large body is present at the radii where the stirred population is observed. In other words, if we consider the stirred population near, say, 150 au, there needs to be at least one large body such that its stirring zone encompasses 150 au. This sets a lower limit on the mass surface density (green upward triangles in Figure 11),

$$\begin{eqnarray}&&{\rm{\Sigma }}\gt \displaystyle \frac{M}{2\pi r{\rm{\Delta }}{r}_{\mathrm{stir}}},\end{eqnarray} \tag{ 14 }$$

where ${\rm{\Delta }}{r}_{\mathrm{stir}}=8\sqrt{3}r{[M/(3{M}_{\star })]}^{\tfrac{1}{3}}$ is the width of the stirring zone (equal to twice the large body's feeding zone, following Equation (3.4) of Ida & Makino 1993). The intersection of the purple line and green triangles represents the maximum mass for a body to stir a given population, and for at least one such body to exist so that its stirred region encompasses a given radius at which the population is observed.

Finally, we can consider the lack of observed radial gaps in the belt as a further constraint on the large body's mass. This is because a body sufficiently large will open a gap of width roughly equal to its chaotic zone (${\rm{\Delta }}{r}_{\mathrm{chaos}}\sim 3r{(M/{M}_{\star })}^{\tfrac{2}{7}}$; e.g., Wisdom 1980; Morrison & Malhotra 2015; Marino et al. 2018). Simply setting the width of the chaotic zone to be less than the resolution of our observations leads to an upper limit to the mass of the bodies within the belt. At 150 au, we find this mass to be <0.15 M_⊕. However, we need to consider the timescale necessary to open this gap. Following Equations (7) and (8) in Morrison & Malhotra (2015), the timescale for a 0.15 M_⊕ body to open a gap at 150 au in the β Pic belt is ∼5.7 Gyr. Then, the lack of radial gaps can be caused not only by the chaotic zone not being resolvable, but also by the large body's gap-opening timescale being longer than the age of the system, where the latter sets a more conservative constraint on the maximum planet mass. At 150 au (50 au), we find that a stirring planet should be less massive than 75 (25) M_⊕ to not produce observable gaps in 23 Myr (brown left-pointing triangles in Figure 11).

In reality, this estimate is uncertain because the clearing timescale of Morrison & Malhotra (2015) is defined as the "time required to reach 50% of the final survivor fraction within the cleared zone" and is therefore not to the timescale to produce an observable gap. The latter may be shorter, for example, if the observing sensitivity is sufficient to detect gap-belt surface brightness contrasts corresponding to a higher survival fraction of particles in the belt. At the same time, the edge-on viewing geometry would make it harder to detect such gap-belt contrast compared to a face-on belt, potentially making the timescale to produce a detectable gap longer. Detailed N-body simulations combined with radiative transfer and observing simulations are needed to accurately pinpoint this timescale (see, e.g., Marino et al. 2018), but are beyond the scope of this work.

Combining the dynamical constraints, we find that the mass of the largest bodies within the belt necessary to stir a hot population such as observed for β Pic over 23 Myr should be between 0.007 and 54 M_⊕ at 150 au. On the other hand, if we consider viscous stirring as the origin of the cold population, we find that the largest bodies must have masses below 0.4 M_⊕ at 150 au (or 0.15 M_⊕ at 50 au) to have maintained such a low dynamical excitation. These masses are all far below and therefore consistent with current detection limits from direct imaging (Lagrange et al. 2018).

To conclude, we underline that viscous stirring by large bodies within the belt alone cannot produce the observed complexity of the vertical structure of the β Pictoris belt. However, it could produce either the hot or cold population, while requiring another mechanism to produce the other. If that were to be the case, the observed dynamical excitation of the hot (cold) populations could be produced through viscous stirring by bodies of mass 0.007–54 (<0.4) M_⊕ at 150 au.

5.2.3. Outward Migration of a β Pic c

A bimodal distribution of inclinations such as observed in the β Pictoris belt was discovered for our own Kuiper Belt (KB) by B01 (after debiasing was taken into account; see also, e.g., Kavelaars et al. 2009; Petit et al. 2011). For both belts, we clarify that a single, broad non-Rayleigh distribution of inclinations cannot be excluded (e.g., Volk & Malhotra 2011); we just follow the original, bimodal functional form of B01 for direct comparison.

Figure 12 shows how the inclination distribution of the β Pic belt (namely, best-fit Model 2) compares with that of the KB from B01. Overall, we find lower inclinations for both the hot and cold inclination populations of the β Pic belt when compared to the KB, indicating an overall lower level of dynamical excitation. The fraction of all particles residing in the cold population are, however, consistent between the β Pic belt and the KB (${20}_{-4}^{+5}$% and ∼19%–26%), with both belts' masses being dominated by the hot population.

Zoom In Zoom Out Reset image size

In the Kuiper Belt, the high inclinations seen for the hot classicals as well as the resonant and scattered populations set important constraints on its formation scenario (for a review, see, e.g., Nesvorný 2015). These populations likely formed in a massive planetesimal disk within 30 au and then got implanted into their present-day orbit by outward-migrating Neptune. The high inclinations can be attained through sweeping resonances (e.g., Levison & Morbidelli 2003) or scattering by Neptune during its migration (Gomes 2003). Nesvorný (2015) showed that the latter mechanism can work as long as Neptune's migration was slow (>10 Myr).

By analogy with the Kuiper Belt, we consider whether β Pic's structure could be produced by outward migration of a hypothetical β Pic c planet located at the inner edge of the belt. The presence of such a planet has been suggested to produce the SW/NE asymmetry observed for the CO gas (Dent et al. 2014; Matrà et al. 2017), as well as for the small dust grains at wavelengths below 24 μm (e.g., Telesco et al. 2005), following models predicting resonant planetesimals to concentrate at specific longitudes relative to the planet (Wyatt 2003). To produce the asymmetry observed for the CO, two conditions need to be met. First, the 2:1(u) resonance trapping probability (where the (u) indicates planetesimals in 2:1 resonance whose libration center increases from 180° during the migration) needs to be sufficiently high, and second, the planetesimal eccentricities need to be excited to a high enough level (we here require e > 0.3; see Figure 6, top row, of Wyatt 2003).

These eccentricities can be excited by resonant forces, causing them to increase as the radial extent of the migration increases (Equation (22) in Wyatt 2003). We find that the ratio between the semimajor axis of planetesimals after and before migration should be ≳1.8 for eccentricities to reach values above 0.3. For example, outward migration of a planet from 31.4 to at least 57.4 au (close to the inner edge of the belt as observed today) could have swept planetesimals in the 2:1(u) resonance from 50 to at least 91 au and produced sufficiently high eccentricities that an asymmetry may be observed. Given this migration must have taken place for at most the 23 Myr age of the system, we estimate a lower limit on the migration rate of 1.13 au Myr⁻¹.

Ensuring the 2:1(u) trapping probability is sufficiently high (we here require it to be above 50%; see Equation (8) and Table 1 of Wyatt 2003) allows us to link this migration rate to the mass of the migrating planet. We find that a planet with mass of at least 28 M_⊕ must have migrated at least ∼25 au to produce an asymmetry as observed for the CO. Therefore, it is possible that outward migration of a super-Neptune-sized planet may have produced the SW/NE asymmetry observed for the CO and for the small dust at short wavelengths. Can this, however, be reconciled with the lack of asymmetry for the mm-sized dust?

Wyatt (2006) showed that the appearance of resonant structure depends on the size of the observed grains, with grains below a critical size D_crit falling out of resonance. This size depends solely on the migrating planet's mass (and some stellar parameters; see Equation (14) in Wyatt 2003). For a planet of at least 28 M_⊕, this size corresponds to 242 μm, although grains up to ∼10 times this size can show significantly smeared resonant structure. The observed emission in the ALMA 1.33 mm data set is dominated by grains of size similar to the observing wavelength or a few times smaller (see, e.g., Figure 6, bottom left in Wyatt 2003). It therefore seems plausible that the resonant structure could be smeared out to the extent that it is undetected in the 1.33 mm ALMA data as well as in the 24 μm (Telesco et al. 2005) and 70–500 μm (Vandenbussche et al. 2010) data sets. Observations at wavelengths even longer than presented here could clarify whether this is the origin of the lack of observed structure at 1.33 mm.

The model predicts the structure of the short-lived, smallest grains and the CO to be asymmetric. This means that the asymmetric resonance structure is only lost, or smeared out, at intermediate grain sizes (corresponding to mid to far-IR/mm wavelengths). Asymmetric resonant structure should therefore be seen for the smallest grains, for the CO, and for the large resonant grains of sizes ≫242 μm. Note, however, that while the distribution of CO and blow-out grains traces the collisional mass-loss rate from the belt (which is proportional to the relative velocity of particles multiplied by the square of their cross section per unit volume), the distribution of the large resonant grains simply traces the particles' cross sections per unit volume. This has the important implication that the SW–NE asymmetry will be much less pronounced and potentially hardly detectable for large resonant grains when compared to the blow-out grains and the CO. This means that the ALMA 1.33 mm and any longer wavelength data could probe resonant grains, but the resonant structure may be too faint for detection.

In addition to the SW/NE asymmetry and the high inclinations observed at the estimated location of the resonant clumps, the outward-migrating planet needs to be able to produce high inclinations all the way out to the outer edge of the disk (i.e., beyond the resonant clumps), as well as the tilt observed for the CO and blow-out grains. A hot population far out in the disk may be produced via outward scattering, as seen for the scattered population of the Kuiper Belt (Brown 2001). At the same time, the inner tilt seen for CO and small grains may be a projection effect caused by two asymmetric, azimuthal resonance clumps in front of and behind the plane of the sky in a non-perfectly edge-on configuration (see, e.g., Matrà et al. 2017, Figure 8, bottom).

Overall, this β Pic c migration scenario could qualitatively reproduce the observed features of the outer belt, but needs to be tested through future planet migration simulations. These should be tailored to reproduce β Pic's gas and dust structure as seen by ALMA, and should take into account the influence that β Pic b would have on such a planet, which we did not consider here.